Setting up the workspace

# ==============================================================================
# * Libraries
# ==============================================================================

library(move)
library(maptools)
library(sp)
library(rgdal)
library(rgeos)
library(lubridate)
library(dplyr)
library(data.table)
library(leaflet)
library(htmlwidgets)
library(adehabitatHR)
library(sf)
library(ggplot2)
library(ggsn)
library(hrbrthemes)
library(DT)
library(tidyr)
library(ggmap)
library(patchwork)
library(cowplot)

source("R/00_fxns.R")

gps_info_out_dir = file.path("results", "tables")
plot_out_dir = "results/maps/gps/eobs"
plot_facet = file.path(plot_out_dir, "facet")

if( !dir.exists(plot_facet)) dir.create(plot_facet, recursive = TRUE)


# ==============================================================================
# * Read in GPS data with flying behaviour + tagging location
# ==============================================================================
# read in the behavior-classified acc and GPS data
acc_gps_behavL = readRDS("data/acc_behavL/yz_acc_gps_behavL.RDS")

gps = lapply(acc_gps_behavL, `[[`, "gps_behav")

gps = lapply(gps, break_nights_fxn)



Find day roosts

# ==============================================================================
# * Find day roosts
# ==============================================================================

# Kilombero roost where tagging took place
kilombero_colony = data.frame(xlim = c(36.98629, 36.99085),
                                    ylim = c(-7.66273, -7.66052))

# Determine day roosts as those gps fixes collected during the day time
get_roosts = function(df, kilombero_colony) {

    is_roost = ifelse(df$DayTime %in% "day" & df$activity %in% "NotFlying",
                      TRUE, FALSE)

    xlim = kilombero_colony$xlim
    ylim = kilombero_colony$ylim
    
    is_kilombero = df$location.long > xlim[1] &
        df$location.long < xlim[2] &
        df$location.lat > ylim[1] &
        df$location.lat < ylim[2]

    new_roost = is_roost & !is_kilombero
    kilombero_roost = is_roost & is_kilombero

    df$roost[ new_roost ] = "new"
    df$roost[ kilombero_roost ] = "kilombero"
    df$roost_date = date(df$timestamp)

    df
}

gps = lapply(gps, get_roosts, kilombero_colony)

# lapply(gps, function(df) table(df$roost))

gps_lines = lapply(gps, make_lines)


if(FALSE) {
    # call and save maps
    htmltools::tagList(
                   lapply(names(gps), map_gps, gps, gps_lines, map_roosts = TRUE)
               )
}


roost_info = lapply(gps, function(df) {
    cat("*** ", df$tagID[1], " ***\n")
    info = df %>%
        filter(roost == "new") %>%
        dplyr::select(tagID, roost, roost_date) %>%
        group_by(tagID, roost_date) %>%
        summarize(n = n()) %>%
        as.data.frame()

}) %>% rbindlist()
## ***  K5309  ***
## ***  K5310  ***
## ***  K5311  ***
## ***  K5312  ***
## ***  K5313  ***
## ***  K5315  ***
## ***  K5317  ***
## ***  K5319  ***
roost_info
##    tagID roost_date  n
## 1: K5309 2016-11-15 30
## 2: K5310 2016-11-17 30
## 3: K5311 2016-11-15 55

Further subset of “not flying” points

# ** Get time differences
gps_tdiff = lapply(gps, function(df) {
    bat_nights = split(df, df$time_cut)
    time_diff = lapply(bat_nights, function(df_night) {
        df_night$time_diff = difftime(lead(df_night$timestamp),
                                      df_night$timestamp, units = "sec")
        df_night$time_diff = as.numeric(df_night$time_diff)
        df_night
    }) 
    time_diff = rbindlist(time_diff)
    as.data.frame(time_diff)
})


# ** Forage points 
# ==============================================================================
    
find_forage_pts =  function(df, min_diff = 30, kilombero_colony) {
    avg_speed = mean(df$ground.speed)
    # avg_speed = mean(df$ground.speed[ df$activity %in% "Flying"])
    sep_m = (avg_speed * min_diff)
    dist = spDists(as.matrix(df[ , c("utm.easting", "utm.northing")]),
                   longlat = FALSE)
    colnames(dist) = 1:ncol(dist)
    rownames(dist) = 1:nrow(dist)
    coloc = data.table(which(dist <= sep_m, arr.ind=T))
    colnames(coloc) =  c("pt1", "pt2")
    coloc_diff = coloc$pt2 - coloc$pt1
    coloc_consec = which(coloc_diff == 1)
    coloc_raw = coloc[ coloc_consec, ]
    coloc_pts = unique(unlist(coloc_raw[ , c(1:2)], use.names = FALSE))
    coloc_time = df$time_diff[ coloc_pts ]
    # coloc_time[ is.na(coloc_time) ] = 9000 # spuriously large number
    coloc_keep = coloc_pts[ coloc_time >= min_diff ]
    df$forage_pts = FALSE
    df$forage_pts[ coloc_keep ] = TRUE
    df$forage_pts[ df$activity %in% "Flying" ] = FALSE
    df$forage_pts[ df$DayTime %in% "day" ] = FALSE

    xlim = kilombero_colony$xlim
    ylim = kilombero_colony$ylim
    
    is_kilombero = df$location.long > xlim[1] &
        df$location.long < xlim[2] &
        df$location.lat > ylim[1] &
        df$location.lat < ylim[2]

    df$forage_pts[ is_kilombero ] = FALSE
    
    # Include distances
    df$dist = NA
    for(i in 1:(nrow(dist) - 1)) {
        df$dist[i + 1] = dist[i, i + 1]
    }
    # Aggregating subsequent foraging points to one group
    x = df$forage_pts
    y = split(x, df$time_cut)
    z = lapply(y, function(v) {
        c(head(v, 1), tail(v, -1) - head(v, -1) == 1)
    })
    z = unlist(z)
    names(z) = NULL
    z = cumsum(c(head(z, 1), tail(z, -1) - head(z, -1) == 1))
    df$forage_cluster = ifelse(x, z, NA)

    forage_time = aggregate(df$time_diff,
                            by = list(df$forage_cluster),
                            FUN = sum, na.rm = TRUE)
    names(forage_time) = c("forage_cluster", "forage_time")
    df = dplyr::left_join(df, forage_time, by = "forage_cluster")
    
}


gps_forage = lapply(gps_tdiff, find_forage_pts, min_diff = 30, kilombero_colony)

# lapply(gps_forage, function(df) table(df$forage_pts))
# lapply(gps_forage, function(df) table(df$forage_pts, df$activity))

forage_info = rbindlist(gps_forage) %>%
    dplyr::select(tagID, forage_pts) %>%
    filter(forage_pts) %>%
    group_by(tagID) %>%
    summarize(n = n()) %>%
    as.data.frame()

forage_info_night = rbindlist(gps_forage) %>%
    dplyr::select(tagID, forage_pts, time_cut) %>%
    filter(forage_pts) %>%
    group_by(tagID, time_cut) %>%
    summarize(n = n()) %>%
    as.data.frame() %>%
    spread(time_cut, n)


if(FALSE){
    htmltools::tagList(
                   lapply(names(gps_forage), map_gps, gps_forage, gps_lines,
                          map_roosts = TRUE, map_forage = TRUE, map_notflying = TRUE)
               )
}

saveRDS(gps_forage, "data/gps_forage.RDS")
gps_forage_df = gps_forage %>%
    rbindlist()

write.csv(gps_forage_df, "data/gps_forage.csv")

Static maps of gps points

# Map gps points
# ------------------------------------------------------------------------------
kilombero_lonlat = data.frame(site = "kilombero",
                           lon = 277974.9,
                           lat = 9152539)

kilombero_wgs = data.frame(site = "kilombero",
                           lon = 36.98726,
                           lat = -7.662083)


kilombero = st_as_sf(kilombero_lonlat, coords = c("lon", "lat"), crs = 32737)

tagIDs = names(gps_forage)

color_palette = c("#7F3C8D", "#3969AC", "#E68310",
                  "#E73F74", "#008695", "#CF1C90",
                  "#f97b72", "#4b4b8f")

tag_colors = color_palette[ 1:length(tagIDs) ]
names(tag_colors) = tagIDs

# ID = tagIDs[1]
small_area_IDs = c("K5313", "K5315", "K5319")
large_area_IDs = "K5310"
medium_area_IDs = tagIDs[ ! tagIDs %in% c(small_area_IDs, large_area_IDs) ]

area_IDs = list(small = small_area_IDs,
                medium = medium_area_IDs,
                large = large_area_IDs)

bats_sf = gps_forage %>%
    rbindlist() %>%
    st_as_sf(coords = c("utm.easting", "utm.northing"), crs = 32737)


map_extents = lapply(area_IDs, function(IDs, bats_sf) {
    bats = bats_sf[ bats_sf$tagID %in% IDs, ]
    bats_sp = as(bats, 'Spatial')
    map_extent = extent(bats_sp) %>%
        as('SpatialPolygons') %>%
        st_as_sf %>%
        st_set_crs(st_crs(bats_sp)) %>%
        st_buffer(3000)
}, bats_sf)
names(map_extents) = c("small", "medium", "large")


basemaps = lapply(map_extents, function(map_extent) {
    extent = map_extent %>%
        st_transform(crs = 4326) %>%
        st_bbox()
    get_map(as.vector(extent),
            zoom=13,
            maptype = 'satellite', source = 'google')
})


p_sm = lapply(small_area_IDs, plot_gps,
              "small", gps_forage, tag_colors, basemaps,
              kilombero_wgs,
              area_IDs, 
              plot_out_dir)

p_med = lapply(medium_area_IDs, plot_gps,
               "medium", gps_forage, tag_colors, basemaps,
              kilombero_wgs,
              area_IDs, 
              plot_out_dir)

p_large = lapply(large_area_IDs, plot_gps,
                 "large", gps_forage, tag_colors, basemaps,
              kilombero_wgs,
              area_IDs, 
              plot_out_dir)

# ** All plots

ind_bat = gps_forage %>% rbindlist

bats_sp = as(bats_sf, 'Spatial')
all_map_extent = extent(bats_sp) %>%
    as('SpatialPolygons') %>%
    st_as_sf %>%
    st_set_crs(st_crs(bats_sp)) %>%
    st_buffer(5000)


all_extent = all_map_extent %>%
    st_transform(crs = 4326) %>%
    st_bbox()

all_basemap = get_map(as.vector(all_extent),
                      zoom=13,
                      maptype = 'satellite', source = 'google')

sc_info = extent(st_transform(bats_sf, crs = 4326))

g = ggmap(all_basemap) +
    geom_point(data = kilombero_wgs, aes(x = lon, y = lat),
                   color = "#08B8FF", size = 4, alpha = 0.65) +
    geom_sf(data = st_transform(map_extents$small, crs = 4326),
            fill = "#6B8278", color = "#6B8278", size = 0.25, alpha = 0.5,
            inherit.aes = FALSE) +
    geom_sf(data = st_transform(map_extents$medium, crs = 4326),
            fill = "#6B8278", color = "#6B8278", size = 0.25, alpha = 0.5,
            inherit.aes = FALSE) +
    geom_path(data = ind_bat, aes(x = location.long, y = location.lat, color = tagID),
              size = 0.75) +
    scale_color_manual(values = tag_colors) +
    theme_ipsum(base_size = 10,
                base_family = "Helvetica",
                plot_title_size = 12,
                grid = FALSE,
                axis = FALSE, ticks = TRUE) +
        theme(
            plot.background = element_rect(fill = '#D1D6AB', colour = '#D1D6AB'),
            legend.position="none",
            plot.margin=grid::unit(c(1,1,1,1), "mm")
    )

    g_scale =

    g +
    ggsn::scalebar(
              x.min = sc_info[1] - 0.001,
              x.max = sc_info[2] + 0.001,
              y.min = sc_info[3] - 0.001,
              y.max = sc_info[4] + 0.001,
              dist = 4,
              dist_unit = "km",
              transform = TRUE,
              model = "WGS84",
              location = "bottomright",
              height = 0.01,
              st.dist = 0.025,
              st.size = 4,
              box.fill = c("white", "black"),
              box.color = "black",
              border.size = 0.25) +
    xlab("Longitude") +
    ylab("Latitude") +
    ggtitle("") 


fn = file.path(plot_out_dir, paste0("all_tags", ".png"))
cowplot::save_plot(fn, g_scale, base_height = 5)

g_scale

Session info

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] cowplot_1.0.0       patchwork_1.0.0     ggmap_3.0.0        
##  [4] tidyr_1.0.2         DT_0.11.5           hrbrthemes_0.6.0   
##  [7] ggsn_0.5.0          ggplot2_3.3.0.9000  sf_0.8-0           
## [10] adehabitatHR_0.4.16 adehabitatLT_0.3.24 CircStats_0.2-6    
## [13] boot_1.3-22         MASS_7.3-51.4       adehabitatMA_0.3.13
## [16] ade4_1.7-13         deldir_0.1-21       htmlwidgets_1.5.1  
## [19] leaflet_2.0.2       data.table_1.12.8   dplyr_0.8.4        
## [22] lubridate_1.7.4     rgeos_0.4-3         maptools_0.9-5     
## [25] move_3.2.0          rgdal_1.4-4         raster_2.9-5       
## [28] sp_1.3-1            geosphere_1.5-10   
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.1         shiny_1.4.0        assertthat_0.2.1   yaml_2.2.1        
##  [5] gdtools_0.2.1      Rttf2pt1_1.3.7     pillar_1.4.3       lattice_0.20-38   
##  [9] glue_1.3.1         extrafontdb_1.0    digest_0.6.23      promises_1.1.0    
## [13] colorspace_1.4-1   plyr_1.8.5         htmltools_0.4.0    httpuv_1.5.2      
## [17] pkgconfig_2.0.3    purrr_0.3.3        xtable_1.8-4       scales_1.1.0      
## [21] jpeg_0.1-8.1       later_1.0.0        tibble_2.1.3       farver_2.0.3      
## [25] withr_2.1.2        magrittr_1.5       crayon_1.3.4       mime_0.9          
## [29] memoise_1.1.0      evaluate_0.14      xml2_1.2.2         foreign_0.8-71    
## [33] class_7.3-15       tools_3.6.0        RgoogleMaps_1.4.3  lifecycle_0.1.0   
## [37] stringr_1.4.0      munsell_0.5.0      compiler_3.6.0     e1071_1.7-2       
## [41] systemfonts_0.1.1  rlang_0.4.4        classInt_0.4-2     units_0.6-3       
## [45] rjson_0.2.20       crosstalk_1.0.0    labeling_0.3       bitops_1.0-6      
## [49] rmarkdown_1.17.1   gtable_0.3.0       codetools_0.2-16   curl_4.3          
## [53] DBI_1.0.0          R6_2.4.1           knitr_1.26         extrafont_0.17    
## [57] fastmap_1.0.1      KernSmooth_2.23-15 stringi_1.4.5      parallel_3.6.0    
## [61] Rcpp_1.0.3         vctrs_0.2.2        png_0.1-7          tidyselect_1.0.0  
## [65] xfun_0.12